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We describe a robust and efficient chain-of-states method for computing Minimum Energy 
Paths (MEPs) associated to barrier-crossing events in poly-atomic systems, which we call the Aeeel- 
eration method. The path is parametrized in terms of a continuous variable t G [0,1] that plays the 
role of time. In contrast to previous chain-of-states algorithms such as the Nudged Elastic Band or 
String methods, where the positions of the states in the chain are taken as variational parameters in 
the search for the MEP, our strategy is to formulate the problem in terms of the second derivatives 
of the coordinates with respect to t, i.e. the state aeeelerations. We show this to result in a very 
simple and efficient method for determining the MEP. We describe the application of the method to 
a series of test cases, including two low-dimensional problems and the Stone-Wales transformation 
in Ceo. 


I. INTRODUCTION 

In the study of molecular, condensed matter or materi¬ 
als systems one is frequently confronted with the need to 
define a transition path for a given atomic rearrangement 
or chemical reaction. This involves specifying a curve in 
configuration space that goes from an initial state of lo¬ 
cal minimum energy, (reactants), to a final one, r^, 
also of local minimum energy (products), that is repre¬ 
sentative of the manifold of actual trajectories through 
which the system could undergo the transition [T]. The 
most obvious and natural way to define such a curve is as 
a Minimum Energy Path (MEP), i.e. a path that fulfills 
the condition of being a minimum of the potential energy 
surface (PES) in the plane perpendicular to the path at 
any point along its length. Equivalently, the MEP is tan¬ 
gent to the PES gradient, and goes through at least one 
saddle point on its way from to r^. 

There are several reasons why the MEP is a useful con¬ 
cept: firstly, as explained above, it gives a clear mathe¬ 
matical definition to the intuitive idea of reaction mech¬ 
anism. The MEP allows to identify the energy barrier(s) 
and possible intermediate states of the transition in ques¬ 
tion. In systems where those barriers are significant (as 
compared to ksT, frequently the case when the transition 
involves the breaking and forming of chemical bonds), 
identifying the relevant MEPs is a pre-requisite to the ap¬ 
plication of Transition-State Theory-based approaches to 
estimate the reaction rate. There are of course situations 
in which the MEP is not such a useful concept. This hap¬ 
pens when there are many competing paths, none of them 
being overwhelmingly dominant [2], which is the typical 
case in soft-condensed matter systems. Path-sampling 
techniques have been developed to estimate transition 
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rate constants specifically in this kind of system [SH- 
Nevertheless, in hard condensed matter and molecular 
systems, the norm is to have transitions that closely fol¬ 
low a well-defined path. This happens e.g. in the diffu¬ 
sion of atoms and defects in solids, either in the bulk or 
at surfaces, in many isomerization reactions in molecules, 
etc. Given the interest in this type of processes, it is 
hardly surprising that many algorithms devoted to find¬ 
ing MEPs have been developed (see [T] and references 
therein), and many practical applications of such algo¬ 
rithms have been reported in the literature. 

There are two strategies that have been frequently em¬ 
ployed in order to identify and locate a MEP. The first 
one starts by locating the first-order stationary point 
(saddle point) that marks the position of the barrier be¬ 
tween the two minima that one wishes to connect through 
the reaction path. This can be done in a number of ways, 
e.g. using a Hessian mode-following algorithm mi, hy¬ 
brid eigenvector following [MI], or the climbing-dimer 
method m- Once the saddle point has been located, 
the MEP can be obtained by following the steepest de¬ 
scent path on either side of the barrier down to the rele¬ 
vant minima [13] . The second strategy, and the one with 
which we will concern ourselves here, attempts to directly 
obtain the full path, usually represented as a string of 
beads or state polymer, in which each bead represents 
a configuration of the entire system displaced along the 
path. Methods of this kind are frequently referred to as 
chain-of-states methods, and some important examples 
are the Nudged Elastic Band (NEB) |TJ [TH [15] and its 
variant, the Doubly-nudged Elastic Band (DNEB) [16], 
the String [TTIITH] and the Ereely-jointed Chain (EJC) [5] 
methods, although there are others (for a review of earlier 
methods of this kind see m- The objective of this family 
of methods is to define a procedure that will cause the 
state polymer to evolve towards the MEP. Not only must 
the converged path fulfill (to within a specified numerical 
accuracy) the conditions for being an MEP; it must also 
retain the states evenly spaced along the chain in order 



2 


to adequately discretize the MEP over its whole length. 
The NEB method achieves this by introducing harmonic 
spring potentials that couple each bead to its two near¬ 
est neighbors along the chain. The configuration of the 
state polymer is then updated by making each bead fol¬ 
low a direction in configuration space that is given by the 
composition of two forces: 

fNEB 

Here = -V£’-L(ri) is the force derived from the PES 
in state i projected onto the hyper-plane perpendicular 
to the path; this term tends to drive the configuration of 
the chain towards the MEP. The second term, is 

obtained from the force due to the harmonic springs, pro¬ 
jected onto the local path tangent. The effect of this term 
is to keep the beads evenly spaced over the length of the 
path. In its original formulation m, the String method 
also uses the first term in Eq. 0 to drive the state poly¬ 
mer towards the MEP; in contrast, this method does not 
use harmonic springs to keep the states evenly spaced, 
but rather uses an interpolation scheme (typically cubic 
spline interpolation [18] ) to parametrize the path and re¬ 
distribute the beads at regular intervals along its length. 
In its more recent, simplified version [18], the full PES 
force is used, as opposed to its path-normal projection, 
to drive the path towards the MEP. Einally, the EJC 
method uses a transformation from Cartesian to hyper- 
spherical coordinates, effectively imposing an even bead 
separation along the chain. Rather than evolving the 
chain in the direction of the normal force along each 
bead, this method minimizes the mis-alignment between 
the force on each bead and the local tangent. 

The NEB and String methods have been very success¬ 
ful, with numerous applications demonstrating their abil¬ 
ity to locate MEPs in complex multi-dimensional sys¬ 
tems. Although there are differences between the two 
(for recent studies comparing them see [IHEO]), they are 
very similar in spirit, with a common denominator being 
the fact that the chain-of-states configuration is evolved 
towards the MEP directly in configuration (coordinate) 
space. This is actually a feature that all chain-of-state 
methods that we are aware of have in common. In this 
work we contend that there is an alternative formula¬ 
tion of the problem in terms of the acceleration vari¬ 
ables, resulting in a very simple algorithm that does not 
require the introduction of spring potentials or otherwise 
re-positioning beads along the chain to ensure an even 
discretization of the path. We term this algorithm the 
Acceleration method. 

The structure of this paper is as follows: in Sec, [ill we 
describe our formalism and strategy for locating MEPs. 
In Sec.[in]we apply the method to a number of test cases, 
namely two simple toy models of reduced dimensionality, 
and a more realistic multi-dimensional problem involving 
an isomerization reaction in Ceo- Einally, in Sec. jlVj we 
review the main features of our method, point out some 
directions for future work and present our main conclu¬ 
sions. 


II. METHODOLOGY 

Our starting point is a parametrization of the path 
between two stationary points on the PES. We will rep¬ 
resent the path as follows: 

r(t) = (1 - t) + u(t). (2) 

Here r(t) is a vector of length d x Nat, with d being 
the space dimensionality (2 or 3 in the examples dis¬ 
cussed in Sec. jlllj ) and Nat the number of atoms in the 
system; t G [0,1] is a reaction parameter, such that 
r(0) = r/^,r(l) = tb, with ta and vb being the given 
start and end configurations, which are stationary points 
(typically minima) of the PES on which we seek to find 
an MEP; and u{t) measures the deviation of the path 
from the linear interpolation and, by construction, must 
fulfill the boundary conditions u(0) = u(l) = 0 . An¬ 
other requirement we impose on u(t) is that its com¬ 
ponents be continuous and twice-differentiable functions 
of t. The objective, then, is to find u(t) for t G [0,1] 
such that Eq. Q is an MEP of the PES between and 
Vb- This will happen when the gradient of the PES at 
any point t, VE[r(t)], is co-linear to the path (whenever 
VE[r(t)] 7^ 0 ), or, in other words, when the gradient 
component perpendicular to the path is zero. Because 
Eq. 0 constitutes an analytical representation of the 
path, for any given trial path we can calculate the path 
tangent, v(t) = dT{t)/dt, i.e. the velocity^ if we view vari¬ 
able t as (fictitious) time. Likewise, we can also calculate 
the acceleration^ a(t) = d?v(t)/dt^. In particular, v(t) is 
important, since it provides us with a criterium for MEP 
convergence (v(t) and VE[r(t)] must be co-linear). As 
we shall see below, a(t) also plays a major role in our 
scheme. Notice that, given the two boundary conditions, 
there is a biunivocal relationship between u{t) and a(t). 

Although Eq. ([^ offers a continuum representation of 
the path, in practical calculations it is necessary to re¬ 
sort to a discrete representation in terms of a set of N 
replicas of the system, r{tn), where tn = nAt, At = 
1/(A — 1), n = 0,1,... A — 1. This amounts to specify¬ 
ing the components of the u vector at N points (including 
the end points) along the path. Discretizing the path in 
this way does not really pose a drawback, as it is always 
possible to effectively recover an analytical representa¬ 
tion by means of an interpolation scheme, or by using a 
set of suitable continuous functions of t to expand the 
components of u. This can always be done provided N 
is not too small. 

Let us now consider the problem of varying the path 
in search of an MEP. As noted in Sec [l| above, previous 
chain-of-states methods use the coordinates of the beads 
along the path as variational parameters in the MEP 
search. As it is well-known [1] , directly optimizing a path 
in terms of bead coordinates will result in a highly wind¬ 
ing path with unevenly distributed beads, and in gen¬ 
eral does not converge towards an MEP. Different strate¬ 
gies can be adopted to avoid this problem (harmonic 
springs coupling neighboring beads in the NEB method. 
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reparametrizing the path at regular intervals, as in the 
String method, etc). In this work, however, we argue 
that the practical difficulties arising from using the coor¬ 
dinates as variational parameters can be very naturally 
overcome by using instead the accelerations, = a(tn), 
as variational degrees of freedom. The idea is simply 
to adjust iteratively in order to drive the path to¬ 
wards a configuration fulfilling the requirement that the 
force perpendicular to the path, = — V^^[r(t^)] = 0 
for n G [0, — 1]. Our method is summarized as follows: 

1. Given an initial configuration of the path (e.g. the 
linear interpolation between start and end configu¬ 
rations, although other choices are possible) and its 
discretization by means of a number N of replicas, 
construct its representation via Eq. § using some 
appropriate interpolation scheme to define the u(t) 
functions (see below). From this representation cal¬ 
culate = v{tn) and for every bead along the 
path. 

2. Calculate the force at each bead position, and 
from it and the path tangent v^, obtain the force 
component perpendicular to the path, i.e. = 

fn - (fn • Vn) Vn, where Vn = ^n/Wnl- 

3. Update the acceleration vector according to: ^ 

an — Af^, where A is a positive numerical parame¬ 
ter, having dimensions of inverse mass, to be suit¬ 
ably adjusted so as to optimize convergence to¬ 
wards the MEP. 

4. By integrating a suitable interpolation a(t) of the 
new an, obtain new vectors Vn and Un- The inte¬ 
gration constants are fixed by the boundary condi¬ 
tions u(0) = u(l) = 0. 

5. Return to step 2 above, and iterate the procedure 
until the path converges to the MEP. 

Before discussing the details of our practical implemen¬ 
tation of the above scheme the following comments are in 
order. Firstly, the need to perform a double integration 
in t to obtain the new path configuration from a^ (step 
4) may be perceived to be a disadvantage of the method. 
However this is not so: the path, and in particular the 
MEP, is generally a smooth, low-curvature trajectory in 
configuration space. It follows that the components of 
u{t) are also smooth, well-behaved and slowly-varying 
functions of t. Therefore, provided N is not too small, 
and a decent interpolation scheme is used, it is possible 
to insure that the integration is performed with suffi¬ 
cient accuracy. A more fundamental reason to work in 
terms of accelerations is discussed at the end of this sec¬ 
tion. Secondly, by viewing the path as a trajectory, and 
t as its time variable, it is easy to see that step 3 above 
changes only the path-normal component of the accel¬ 
eration. This component affects only the shape of the 
path, i.e. the direction of its tangent vector, v(t), but 
not its modulus. It follows that images are not caused 


to slide up or down the path in any significant way, and 
thus the inter-bead spacing will (to first order) remain 
even. Nevertheless, inter-bead spacing will become un¬ 
even over a sufficiently large number of iterations of the 
scheme due to curvature effects. If needed, the tangential 
components of the acceleration can be scaled by a factor 
smaller than 1, so as to gradually reduce their value dur¬ 
ing the iterative process, which will ensure even spacing 
of the beads. In the examples that follow we found this to 
be unnecessary, although we did it in the second example 
for illustrative purposes. Thirdly, the parameter A intro¬ 
duced in step 3 above determines the rate of convergence 
of the method, and choosing it well is therefore impor¬ 
tant. In the illustrative examples discussed in Sec. |m] 
we have for simplicity adopted the strategy of taking it 
as a constant value, giving overall adequate convergence. 
However better convergence rates may be achieved by al¬ 
lowing A to vary and using e.g. a Hessian-up date scheme 
to choose A appropriately at each step and/or for each 
bead independently. This issue will be the subject of 
future research. 

The general scheme described above could be imple¬ 
mented in a number of different ways; all that is needed 
is a flexible interpolation scheme that allows to construct 
a representation of the u(t) vector components from the 
bead positions, or rather, their second derivative with 
respect to t (which enter in the acceleration), and to per¬ 
form the reverse process of integrating 3.new{t) to ob¬ 
tain the new configuration of the path (step 4). This 
could be done e.g. using cubic spline interpolation, or 
any other suitable interpolation scheme. In particular, 
we have found a Fourier sine series representation of the 
u vector components to be particularly convenient. In 
our implementation we represent them as follows: 

N-2 

u(^) = E UnSin(w„t), (3) 

n=l 

where = nir. The N — 2 nonzero Fourier coefficients 
Un are fixed by the N — 2 nonzero values u^, with 0 < 
n < N — 1. Eq. ^ obeys the boundary conditions u(0) = 
u(l) = 0 by construction. Another advantage is that the 
first and second derivatives u'(t), and u"(t), are similarly 
given as cosine and sine series, respectively. It is therefore 
very simple to obtain u{t) and v(t) from a(t), as required 
by step 4 of our algorithm. Indeed, following step 3 one 
obtains new accelerations a(t) = u"(t), which, by virtue 
of Eq. have components of the form 

N-2 

a{t) = y; a„sin(w„i), (4) 

n=l 

where again the Fourier coefficients a^ can be obtained 
from a.n by Fourier transform techniques. Now, Eq. 0 
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can be integrated two times to give 


JV-2 - 

W —COs{uJnt) + Co, 
n—1 

(5) 

N-2 ^ 

W ^sin(w„i) + Cl. 

(6) 


The boundary conditions fix the values of the integration 
constants to be Cq = and Ci = 0. 

In the next section we will show that the method just 
described is robust, efficient and stable. Before we de¬ 
scribe its application to specific examples, it is worth 
pausing to reflect on the reasons for its stability. One 
may naively assume that a similar scheme to ours, but 
formulated in terms of coordinates instead of accelera¬ 
tions (i.e. using ^ -1- Af^ in step 3) should work 

just as well, thus obviating the need to integrate acceler¬ 
ations to obtain the coordinates and velocities. However 
practical experience shows that this is not the case, as is 
well documented [T] . Such a scheme results in a snake-like 
path dominated by high-frequency error components that 
never converges to the MEP. This, however, does not hap¬ 
pen in our scheme, and in Eq. we can see the reason 
for this: in the integration step to obtain the coordinates 
each acceleration component is scaled by the inverse of 
its corresponding frequency squared, thus effectively act¬ 
ing as a filter to high frequency error components. As 
a consequence, the path evolves more smoothly towards 
the MEP, which allows a faster convergence without de¬ 
veloping kinks or twists in the process. 


III. RESULTS 

In order to demonstrate the effectiveness of the 
methodology presented in Sec. [H] above, we describe here 
its performance in three specific cases of MEP location. 
The first two examples we consider are simple 2D poten¬ 
tial energy surfaces (PES), namely the modified LEPS 
potential (model II in [T]) and the Muller potential [22]. 
As a third example we consider the multi-dimensional 
problem of the Stone-Wales isomerization transition [23] 
between the (buckminsterfullerene) and the C 2 v low¬ 
est energy isomers of Ceo- The first two examples are 
simple toy models of reduced dimensionality, but never¬ 
theless contain all the essential ingredients of the problem 
in the more general, multi-dimensional case. In spite of 
their simplicity and 2D character, they constitute chal¬ 
lenging test cases for any methodology that aims to be a 
viable alternative for the location of MEPs. In all cases 
discussed below we took as initial guess for the MEP a 
simple linear interpolation between the end points of the 
path, which invariably were chosen as two previously lo¬ 
cated minima on the corresponding PES. The number 
of beads or replicas of the system along the path was 
varied between a minimum of 10 and a maximum of 30, 




FIG. 1: (Color online) The top panel shows the MEP obtained 
for the LEPS potential using 30 replicas of the system along 
the path. The black cross marks the position of the saddle 
point between the two minima, and the minima, marking the 
start and end of the path, are labeled as A and B. For com¬ 
parison, the lower panel displays simultaneously three MEPs 
obtained with 10, 20 and 30 beads. 

although individual tests have been also made with bead 
numbers outside this range. 


A. Modified LEPS potential 

The modified LEPS potential (see [T] for details of its 
definition) possesses two minima, the first one of which 
is located at x = 0.7415,^ = 1.3034, will be labeled 
as A in what follows, and is the global minimum on 
this PES. The second, local minimum B, is found at 
X = 3.0012,^ = —1.3040. A barrier separates the val¬ 
leys of each minimum, with a saddle point located at 
X = 2.021,^ = —0.173. A contour plot of this potential 
is shown in Eig. 0 . Also shown in the figure is the con¬ 
verged MEP [21] that resulted with 30 beads in the path 
discretization. As expected, the obtained MEP cuts per¬ 
pendicularly the PES contour lines, and passes through 
the saddle point at the top of the barrier between the 
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FIG. 2: (Color online) The energy profile of the LEPS po¬ 
tential along the MEP obtained with 10, 20 and 30 replicas 
of the system along the path, t = 0 corresponds to the first 
minimum [A, see Eig. 01. and t = 1 to the second one (B). 

A and B valleys (marked with a black cross on the fig¬ 
ure). Note how the replicas (shown as red dots) remain 
roughly evenly spaced along the MEP, in spite of the fact 
that no harmonic springs are used to impose this, in con¬ 
trast to the case of the NEB method; neither have they 
been artificially redistributed, as in the String method. 
Even the top of the barrier remains well described by a 
sufficient density of beads. 

It is worth noticing that in spite of the simplicity and 
reduced dimensionality of this PES, the MEP has sharp 
bends, where the path turns by nearly 90 degrees (when 
climbing out of the A valley and down into the B val¬ 
ley). Such regions of high curvature would pose a chal¬ 
lenge for any simplistic approach to MEP location, but 
our methodology encounters no particular difficulty with 
these regions. 

Fig. (Il displays the energy profiles along the MEP 
when the latter is discretized with 10, 20 and 30 beads. 
As can be seen there, using only 10 beads results in a 
relatively rugged description of the MEP, although the 
general features of the path, such as the barrier height, 
are reasonably well reproduced even in this case. With 
20 and 30 beads a much smoother and accurate repre¬ 
sentation of the MEP is obtained, as evidenced from the 
fact that the energy profiles in these cases are hardly dis¬ 
tinguishable on the plot. The description of the barrier 
summit is slightly more accurate with 30 replicas due 
to the increased density of beads, but elsewhere the two 
profiles are practically identical. 

B. Muller potential 

Let us now consider the case of the Muller poten¬ 
tial [22]. In contrast to the LEPS model seen above, this 
PES has three minima, and two saddle points separating 




EIG. 3: (Golor online) The MEP obtained for the Muller po¬ 
tential using 30 replicas of the system along the path. The 
black crosses mark the position of an intermediate local min¬ 
imum and of two saddle points between the two end minima. 
The minima are labeled as A (the absolute minimum), B, 
and the intermediate one, G. Two saddle points mark the 
position of the barriers separating the minima, the first one 
being labeled as Ae^G, and the second one labeled as Ge^B. 
The upper panel shows the MEP obtained without perform¬ 
ing any scaling of the tangential acceleration components; as 
can be seen, the bead distribution becomes somewhat uneven. 
In the lower panel, scaling the tangential acceleration compo¬ 
nents by a factor of 0.99 at every convergence iteration results 
in an even distribution of beads along the MEP. 


them. Although still only a 2D model, the presence of 
more stationary points on the PES constitutes an added 
challenge for MEP location algorithms. A contour plot 
of this PES is shown in Fig. (§, together with the lo¬ 
cation of the various stationary points. The minima are 
labeled as A (the absolute minimum), with coordinates 
X = —0.558,^ = 1.442, B, x = 0.623,^ = 0.028, and the 
intermediate one, C, located at x = —0.05,^ = 0.467. 
Two saddle points separate the valleys corresponding to 
each minimum, at x = —0.822,^ = 0.624 (labeled as 
Ae^C) and x = 0.212,^ = 0.293 (labeled as Ce^B), re¬ 
spectively. Like in the case of the LEPS model, we con¬ 
sidered path discretizations using 10, 20 and 30 beads. 
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FIG. 5: (Color online) Scheme of the Stone-Wales transition. 
The central bond, highlighted in red, rotates by 90°; for this 
to happen two bonds have to be broken, such as the ones 
marked by dashed lines, and re-formed after the rotation of 
the central bond takes place. 


FIG. 4: (Color online) The energy prohle of the Muller po¬ 
tential along the MEP obtained with 10, 20 and 30 replicas 
of the system along the path, t = 0 corresponds to the first 
minimum [A, see Fig. §1. and t = 1 to the second one (B). 

and in each case the initial path was taken as the linear 
interpolation between A and B. In the upper panel of 
Fig. (§ we plot the resulting MEP m using 30 repli¬ 
cas, without using any scaling of the tangential acceler¬ 
ation components to keep the beads evenly spaced. As 
expected, the obtained MEP goes through the intermedi¬ 
ate minimum (C) and the two saddle points located along 
the path. While the finest description of the MEP is ob¬ 
tained with 30 replicas, coarser path descriptions with 
10 and 20 beads (not shown) also track the MEP cor¬ 
rectly. The converged MEPs have beads roughly evenly 
distributed along the whole length, even without scaling 
the tangential components of the acceleration. For com¬ 
parison, in the lower panel of Fig. we show the MEP 
obtained when the tangential acceleration components 
are scaled by a factor 0.99 at each step of the iterative 
process; as can be seen, the MEP that results in this case 
has homogeneously distributed beads. 

Fig. @ shows the energy profile along the converged 
MEPs with the different path discretization used in the 
calculation. As can be seen by comparing Figs. Q 
and the energy profile on the Muller PES has more 
features (two peaks and an intermediate valley) than that 
of the LEPS model. Using onlylO beads to describe the 
MEP results in a rather coarse description of the path, 
but even at this level all features of the energy profile are 
captured and reasonably described. 


C. Stone-Wales transformation in Ceo 

As a final example we will consider the case of the 
Stone-Wales (SW) transformation [23] in Ceo- It is 
known that Ceo has 1812 fullerene isomers (i.e. cage 
structures formed exclusively by 20 hexagons and 12 pen¬ 
tagons). Only one of these isomers obeys the isolated 


pentagon rule [24]; it is the Ih structure known as buck- 
minsterfullerene, having every pentagon surrounded by 
five hexagons (there are no pentagon-pentagon adjacen¬ 
cies in this structure). This is the most stable structure 
of Ceo; in all other Ceo isomers pentagon-pentagon ad¬ 
jacencies are present, incurring an energy penalty that 
renders these structures less stable than the Ih isomer. 
Stone and Wales [23] were the first to point out that it 
was possible to transform a given fullerene isomer into 
a different one by means of the rotation of a C-C bond 
connecting two hexagons and two pentagons. The rota¬ 
tion of such a bond around its center interchanges the 
positions of the hexagons and pentagons, as illustrated 
in Fig. <§• Given the importance of the SW transforma¬ 
tion in the growth of carbon nanostructures m and as a 
stress-release mechanism in carbon nanotubes [ 2 a[ 27 ],it 
has been extensively studied at different levels of theory 
(see e.g. ref. [28l[29] and references therein). 

Here we will consider the SW transformation between 
the Ih buckminsterfullerene and the C 2 v lowest-energy 
isomers of Ceo- The energetics of the system has been 
described with three different models, namely the many- 
body potential due to Tersoff m, the orthogonal tight- 
binding (TB) model of Xu et al [31], known as MDTB, 
and the non-orthogonal TB model model due to Porezag 
et al. [32], known as DFTB. All these models have been 
extensively used in the study of carbon-based systems. 
To perform the calculations described below we have cou¬ 
pled the MEP-search method described in Sec. |TI| with 
the Trocadero code [33] , which contains implementations 
of all these energy models. The connecting path between 
the Ih and C 2 v isomers has been discretized with 30 repli¬ 
cas; tests with different numbers of replicas (both smaller 
and larger than 30) were also carried out, leading to prac¬ 
tically indistinguishable MEPs. In Fig. we plot the 
energy profiles along the MEP that resulted from these 
calculations. The energy of the ground state (isomer Ih) 
is taken as the zero of energy. Let us first discuss the 
MEP that results with the Tersoff potential. According 
to this model, the energy difference between the Ih and 
C 2 v isomers is only 0.75 eV, much below typical values 
found with first-principles electronic structure methods. 
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FIG. 6: (Color online) The energy profile along the Stone- 
Wales transition MEPs in Ceo- The structures shown cor¬ 
respond to the Ih isomer (buckminsterfullerene, the ground 
state) on the left, the C 2 v isomer on the right, and the sad¬ 
dle point conhguration between the two, as calculated with 
the DFTB model. MEPs obtained with the Tersoff potential 
and the MDTB and DFTB tight binding models are shown 
for comparison. Red and blue arrows mark the energy of the 
DFTB and MDTB saddle points, respectively (see text). 


which are about twice as large [28l [29]. Furthermore, 
the MEP obtained from the Tersoff description is very 
unusual, showing two sharp peaks with discontinuous 
derivative either side of a local minimum at the center 
of the barrier. The peaks do not become smoother if 
a larger number of beads is used to discretize the path; 
they appear to be features of this model potential. The 
local minimum found close to the center of the barrier 
appears because this potential tends to over-stabilize the 
carbon atoms at the edge of the SW motif in the central 
configuration displayed at the top of Fig. turning 
the structure into a local minimum, instead of the saddle 
point that one would expect to find close to the barrier 
center. The maximum barrier height found along the 
Tersoff MEP is 5.44 eV, somewhat lower than that found 
by Marcos et al. [34] (5.58 eV) using the same model; 
these authors reported a low-energy path from the Ih to 
C 2 v isomers by linear interpolation between a series of 
intermediate minima that they identified from molecular 
dynamics simulations. The overall shape of their energy 
profile is similar to the MEP reported here, but contains 
actually three local minima, not counting the end struc¬ 
tures. We believe that the two local minima on either 
side of the barrier they report are actually not there, and 
are seen in their energy profile because their path is not 


a true MEP. Marcos et al. [34] attached some physical 
significance to the local minima they found, arguing that 
they would facilitate the global SW transition, as each 
intermediate step is subject to lower barriers than the 
overall process. We rather believe that the local mini¬ 
mum at the center of the barrier (the only one we find) 
is an artifact that results from the poor transferability of 
the Tersoff potential to situations far departed from the 
configurations considered in its parameter fitting. 

The main difference between the Tersoff potential and 
the TB models discussed next is that the latter incorpo¬ 
rate a description of the (valence) electronic structure of 
the carbon atoms, albeit at a semi-empirical level [35] . 
The TB models themselves are similar in spirit, but dif¬ 
fer mainly in the fact that the MDTB model assumes 
the underlying basis set to be orthonormal, while the 
DFTB model explicitly incorporates their overlap, which 
in principle makes the model more transferable. These 
TB models provide a much more credible picture of the 
SW transition, that is in better agreement with what is 
predicted by higher levels of theory. Indeed, the MDTB 
model predicts the energy difference between the end iso¬ 
mers to be 1.41 eV, while the DFTB model gives a value 
of 1.7 eV, both much closer to the range of values result¬ 
ing from first-principles calculations, which average at 
about 1.8 eV [28l|29]. Both models also provide similar 
barrier heights: 6.45 eV (MDTB) and 6.44 eV (DFTB), 
which again are close to values predicted by higher levels 
of theory, averaging at about 7.5 eV. As can be seen in 
Fig. (§, none of the TB models predicts the existence 
of any local minimum along the MEP, in contrast to the 
Tersoff potential. The obtained paths are smooth, with¬ 
out sharp features. In spite of the higher dimensional¬ 
ity of this problem {ZN — 6 = 174 degrees of freedom), 
compared to the LEPS and Muller potentials considered 
above, it is still the case that the replicas remain roughly 
evenly spaced along the converged MEPs, even in the 
case of the Tersoff potential description, in the presence 
of the sharp features observed there. 

The red and blue arrows on either side of the bar¬ 
riers shown in the main panel in Fig. §) mark the en¬ 
ergy of the transition-state configuration according to the 
DFTB and MDTB models, respectively. These configu¬ 
rations have been located using the mode-following Ra¬ 
tional Function Optimization (RFO) approach of Baner- 
jee et al. laiH]. This algorithm searches for stationary 
points of the PES using both first-derivative (gradient) 
and of second-derivative (Hessian) information. A saddle 
point can be located following a chosen eigenvector of the 
Hessian up-hill in energy, while minimizing along all the 
other modes. In our case we maximized along the eigen¬ 
vector associated to the lowest Hessian eigenvalue, start¬ 
ing from the highest energy configuration found along 
the MEP. The search for the saddle point was assumed 
to have converged once the resulting structure had a 
gradient with components smaller than 10“^ eV/A and 
the Hessian had a single negative eigenvalue. The final 
DFTB configuration is illustrated as the central structure 







at the top of Fig. (^; the one obtained with the MDTB 
model is very similar and is not shown. As can be seen 
in the figure, the energies of the transition states located 
in this way are very close to the barrier heights predicted 
from the MEPs, indicating that these are well converged. 

In spite of the fact that the DFTB and MDTB mod¬ 
els predict comparable C 2 v — energy differences and 
barrier heights, their MEPs have differences, as well as 
similarities. The MDTB path results in a much narrower 
barrier. Analyzing the structures along the MDTB path, 
one can see that as the path moves away from the 
start configuration, the Ceo sphere distorts to become an 
ovoid before any SW rotation happens; the sharp energy 
rise that results in the barrier occurs only once chemical 
bonds break to allow for the dimer rotation to occur. In 
the case of the DFTB model, however, dimer rotation 
begins much sooner, resulting in a wider barrier. This 
is consistent with the fact that the MDTB model pre¬ 
dicts much softer vibrational frequencies than the DFTB 
model for Ceo- 

The saddle point configuration closest to the maximum 
of the MEP shown in Fig. (§ has C 2 point-group sym¬ 
metry. It is very similar to the structure with the same 
symmetry reported by Bettinger et al [29] on the ba¬ 
sis of DFT calculations. The same authors reported a 
second, asymmetric, transition state involving a carbene 
intermediate. We do not find such a structure with either 
TB model used in this work. It is interesting that Walsh 
and Wales [36] found a different saddle point configura¬ 
tion to the one we obtain, even though in principle they 
used the DFTB model. Their structure is asymmetric, 
with the rotating C-C bond highlighted in red in Fig. ^ 
tilted towards one side of the open cage, forming a trian¬ 
gle with the under-coordinated atom at the vertex of the 
nearby hexagon (see Fig. 1 of their paper). We have also 
searched for this structure using the RFO miH] method, 
but have not been able to locate it. Starting the RFO 
search from a structure similar to theirs converges to the 
same symmetric C 2 structure that we find from our cal¬ 
culated MEP. The reason for this discrepancy remains 
unclear, but it is most likely due to the use of different 
parametrizations of the same TB model in their work and 
ours. 


IV. CONCLUSIONS 

We have presented a new algorithm for the search of 
minimum energy paths (MEPs) between two given sta¬ 


tionary points on a potential energy surface, and demon¬ 
strated its use and robustness when applied to a number 
of low-dimensional toy models and to the Stone-Wales 
transformation in Cqq. Our method falls in the class of 
algorithms known as chain-of-states methods^ requiring 
only gradient information from the PES. But in con¬ 
trast to well-known examples of other methods of this 
kind, such as the Nudged-Elastic Band^ it does not rely 
on the introduction of additional force terms acting on 
the states, nor does it rely on redistributing the states 
along the path, as in the String method, to converge to 
the MEP. In spite of this, in all the test cases we have 
considered the method converges smoothly, retaining the 
images discretizing the path approximately evenly spaced 
over its length. Because the method uses an analytical 
representation of the path, it is possible, if desired, to 
refine the path discretization by inserting new images 
where required, or to scale the tangential acceleration 
components to ensure an even spacing, although we stress 
that in the test calculations we have reported above this 
was found to be unnecessary. The method has an appeal¬ 
ing simplicity, making it easy to implement and combine 
with existing atomistic simulation codes. 

Our practical implementation of the scheme is robust 
and efficient. Nevertheless it is susceptible to improve¬ 
ment in a number of ways. In particular, the acceler¬ 
ation update scheme (step 3 of the method) is akin to 
a steepest-descent method. Previous experience in the 
NEB and String methods has shown that more sophisti¬ 
cated update schemes can significantly improve the rate 
of convergence towards the MEP, and we expect this to 
be the case here as well. We will explore this issue in 
future research. We also hope to demonstrate the useful¬ 
ness of the present scheme in a wider set of case studies. 
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